Elastic wave propagation in confined granular systems 
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We present numerical simulations of acoustic wave propagation in confined granular systems 
consisting of particles interacting with the three-dimensional Hertz-Mindlin force law. The response 
to a short mechanical excitation on one side of the system is found to be a propagating coherent 
wavefront followed by random oscillations made of multiply scattered waves. We find that the 
coherent wavefront is insensitive to details of the packing: force chains do not play an important role 
in determining this wavefront. The coherent wave propagates linearly in time, and its amplitude 
and width depend as a power law on distance, while its velocity is roughly compatible with the 
predictions of macroscopic elasticity. As there is at present no theory for the broadening and decay 
of the coherent wave, we numerically and analytically study pulse-propagation in a one-dimensional 
chain of identical elastic balls. The results for the broadening and decay exponents of this system 
differ significantly from those of the random packings. In all our simulations, the speed of the 
coherent wavefront scales with pressure as p^^^; we compare this result with experimental data on 
various granular systems where deviations from the p^^^ behavior are seen. We briefiy discuss the 
eigenmodes of the system and effects of damping are investigated as well. 

PACS numbers: 45.70.-n, 43.40.+S, 46.40.Cd, 46.65. -fg 



I. INTRODUCTION 

The behavior of granular systems is often strongly in- 
fluenced by fluctuations and inhomogeneities on the scale 
of the granular particles, which precludes or complicates 
the derivation of macroscopic laws from grain-level char- 
acteristics. In quasi-static granular packings the strong 
fluctuations in the intergrain forces are a striking exam- 
ple of this heterogeneity The fact that the grain- 
grain contacts which support large forces are usually cor- 
related in a line-like fashion over distances of several par- 
ticle diameters leads to so called force chains. (In static 
equilibrium a contact transmitting a large force is of- 
ten balanced with a single contact on the opposite side 
of the grain, and this is repeated on several subsequent 
grains.) These are clearly visible in experiments on two- 
dimensional (2D) packings using photoelastic disks IJ 
and in simulations, even though their precise definition 
is not agreed on at present. 

More importantly, it is not clear what properties of 
the granular media are affected by these force chains or 
indeed by the broad distribution of interparticle forces. 
It is under debate whether granular media can be de- 
scribed as elastic, and since a static granular packing is 
a quenched system, an important issue is at what length 
scale such a continuum (elastic) theory would become 
appropriate for a granular medium. 

An important quantity that can probe these issues is 
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the propagation of sound waves. Surprisingly, even the 
scaling with pressure of a continuum quantity like the 
speed of sound is still a matter of debate — even though 
clearly the average interparticle force scales with the pres- 
sure, the number of interparticle contacts also increases 
with pressure, and it is being debated whether this might 
affect the dependence of the speed of sound on pressure 
QtHil^' Moreover, one might argue as follows for an in- 
terplay between force chain type correlations and sound 
propagation. In a simple ID array of grains the group 
velocity of waves is higher for a more compressed array 
because of the nonlinearity of the force law Q ■ Similarly, 
a simple calculation shows that in a continuum elastic 
medium with a hard block embedded in soft surrounding, 
the small amplitude acoustic waves propagate principally 
in the hard material. From these observations one might 
wonder whether the acoustical waves that travel through 
the granular medium first are transmitted mostly by the 
strong force chains. If so, one might be able to extract de- 
tailed information about the force chains from the trans- 
mitted acoustic signal. 

In this paper we explore numerically and analytically 
what ultrasound experiments Q might tell about the mi- 
croscopic and mesoscopic features of granular media: can 
they tell us anything about the existence of force chains? 
Can ultrasound probe on what length scale we can view 
a granular medium as an almost homogeneous random 
medium? Are there generic features of the propagation 
of a short pulse which we can uniquely relate to the ran- 
dom nature of a granular pack? These are some of the 
issues we have in mind in this paper. 

Using numerical simulations in 2D and 3D, we show 
that sound does not predominantly travel along force 
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chains: the leading part of the wave, which propagates 
after being excited by a short pulse at one end of the 
medium, is better characterized as a propagating rough 
front, like what was found in bond-diluted models [3- 
Similarly to experiments of this type, we can separate 
the transmitted acoustic signal into an initial coherent 
part and a subsequent random part. The coherent wave 
propagates essentially linearly in time, defining a time- 
of-flight sound velocity. We have calculated the effective 
elastic constants for our packings, which under the as- 
sumption that the granular packing can be viewed as a 
continuum directly lead to a simple expression for this ve- 
locity. Surprisingly, the observed time-of-flight velocity 
of the coherent wave is roughly 40% larger than predicted 
by continuum elasticity. 

We also study the scaling of sound velocity and elastic 
constants with pressure. For a fixed contact network with 
Hertzian forces which stay proportional to the pressure p, 
the elastic constants scale as p^^^ and the sound velocity 
as p^/^ . In the pressure range studied here we find no dis- 
cernible deviation from the p^^^ behavior for the time of 
flight velocity; the situation for the velocities calculated 
from the effective elastic constants is more convoluted. 

Discrete numerical simulations have been used to eval- 
uate the elastic moduli of packings of spherical beads 
before 0, H, ) from which continuum sound speeds can 
be deduced. However, to our knowledge, wave propaga- 
tion has never been directly addressed numerically in a 
granular material. Measurements of overall elastic prop- 
erties do not probe the material on the mesoscopic scale 
and overlook potentially interesting properties such as 
dispersion and attenuation. Hence the interest and mo- 
tivation of the present work, which directly copes with 
the properties of a traveling ultrasonic pulse in a model 
granular packing. 

We find that there are various other interesting aspects 
of the problem of pulse propagation in granular systems 
which appear to have received little attention so far: for 
disordered systems the amplitude of the coherent wave 
decays as a power law as it propagates, while its width 
increases linearly. As we are unaware of systematic stud- 
ies of these issues so far, we also consider, as a reference, 
the propagation of a sound pulse in ID chains and in 2D 
triangular lattices of identical elastic balls. We show that 
also here both the amplitude and width of the coherent 
wave behave as power law of the distance. We calculate 
both these exponents and the wave front analytically, and 
show that the broadening of the pulse and the decay of 
the width is much slower than in the 2D disordered sys- 
tem. A more detailed experimental and theoretical study 
of these aspects might therefore yield an important way 
to probe granular media in the future. 

This paper is organized as follows. First we review 
in Sec.nisome of the relevant results of experiments and 
numerical simulations. Then Sec. lHII describes the details 
of our numerical model: how the packings are obtained 
and how the small amplitude oscillations are analyzed. In 
Sec. II VI we present qualitative and quantitative results of 



our simulations. In Sec. E] we compare these results to 
theoretical models of a ID chain and a triangular lattice 
of identical balls. Finally Sec. I VII concludes the paper. 



II. ELASTICITY AND SOUND PROPAGATION: 
EXPERIMENTAL RESULTS 

There have been a number of related experiments that 
considered sound propagation in granular systems. Be- 
cause of their direct relevance we briefly review their 
main results. It is useful to keep in mind that, in prin- 
ciple, two types of experiments can be performed: either 
one drives the system with constant frequencies and fo- 
cuses on spectral properties, or one drives the system 
with short pulses, testing propagative features. Since 
wave propagation is traditionally described in the frame- 
work of macroscopic linear elasticity, we also briefly evoke 
some of the measurements of elastic moduli as obtained 
in laboratory experiments. 

Liu and Nagel 0,0, 0| studied acoustic sound prop- 
agation through an open 3D granular assembly. They 
prepared a 15-30 grain diameter deep layer of glass beads 
in an isolated box. In their setup the top surface was 
free and subject to gravity only. The sound source was 
a vertical extended plate, embedded within the granu- 
lar layer, and the detectors were accelerometers of size 
comparable to the grains. They identifled three distinct 
sound velocities. From the response to a short pulse, 
the ratio of the source-detector distance and the time of 
flight gives Ctof = -^/Tflight = 280 ± 30 m/s, while the 
dependence of the time of maximum amplitude of the 
response on source-detector distance yields Cmax.rosp = 
dL/dTniax — 110 ± 15 m/s. In case of harmonic excita- 
tion, the group velocity deflned by the frequency depen- 
dence of the phase delay is Cgroup = ^TTLdv/dcf) = 60 ± 10 
m/s. From these incompatible values they concluded 
that the granular packing cannot be considered as an 
effective medium for sound propagation, as the trans- 
mission is dominated by the strong spatial fluctuations 
of force networks. A consequence of this is the extreme 
sensitivity on small changes (e.g., heating one bead by 
less than a degree 0), and the power spectrum on 
response to harmonic excitation. Many aspects of these 
experiments have been conflrmed by numerical simula- 
tions in a model system of square lattice of non-identical 
springs il2(| and with Hertzian spheres 0| . 

Jia and coworkers measured sound propagation 

through a conflned 3D granular system. They flUed a 
cylinder of radius and length 15-30 grain diameters with 
glass beads, compacted by horizontal shaking, closed off 
with a piston, and applied an external pressure on the 
piston. The sound source was a large flat piezo trans- 
ducer at the bottom of the cylinder, and the (variable 
size) detector was at the top wall. They measured the 
response function for a short pulse, and observed that it 
can be divided into an initial coherent part, insensitive to 
the details of the packing, followed by a noisy part, which 
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changed significantly from packing to packing. The ratio 
of the amphtudes of the two parts of the signal depended 
on the relative size of the detector and the grains, and 
also on additional damping (e.g., wet grains 

Gilles and Coste jlTj studied sound propagation in a 
confined 2D lattice of steel or nylon beads. They ar- 
ranged the beads on a triangular lattice of hexagonal 
shape, 30 beads on each side, and isotropically com- 
pressed the system. This resulted in a regular array of 
bead centers, but irregular intergrain contacts because of 
the small polydispersity of the beads. Both the sound 
source and the sensor were in contact with a single bead 
at opposite sides of the hexagon. They also found that 
the response to a short pulse was an initial coherent sig- 
nal, followed by an incoherent part. While the whole 
signal was reproducible for a fixed setup, they observed 
a high correlation factor between the coherent signals of 
different packings, and low correlation between the inco- 
herent parts. 

Apart from these recent experiments carried out in the 
physics community, it is worth recalling that many inter- 
esting measurements of the elastic or acoustical charac- 
teristics of granular materials are to be found in the soil 
mechanics and geotechnical engineering literature ( ^3 
is a recent review stressing the need for sophisticated 
rheological measurements of soils, including elastic prop- 
erties) . The elastic properties of granular soils have been 
investigated from quasistatic stress-strain dependencies, 
as measured with a triaxial or a hollow c ylin der po| 
apparatus, by "resonant column" devices |2ll|22| (which 
measure the frequency of the long wavelength eigenmodes 
of cylindrically shaped samples), and from sound prop- 
agation velocities IS, ,^0^ 2^ -23. l25l| . One remarkable 
result is the consistency of moduli values obtained with 
various techniques, provided the applied strain incre- 
ments are small enough (i.e., typically, lower than 10~^). 
Thus, the agreement between sound propagation and res- 
onant column results is checked e.g. in |2^, while po| 
shows consistency between sound propagation velocities 
and static moduli for very small strain intervals. 

In that community, wave velocities are most often de- 
duced from signal time-of-flight measurements between 
pairs of specifically designed, commercially available 
piezoelectric transducers known as bender elements [23i . 
l2g | which couple to transverse modes or bender- extender 
elements, which also excite longitudinal waves p^ . The 
typical size of such devices is 1cm, sand grains di- 
ameters are predominantly in the 100/im range, a typi- 
cal propagation length (specimen height or diameter) is 
10cm, and confining stresses range from 50 or 100 kPa to 
several MPa. Those experiments are thus comparable to 
the ones by Jia the material being however probed 
on a somewhat larger length scale. The shape o f sig nals 
recorded by the receiver are similar (see e.g. p3 l25j'). 

The soil mechanics literature on wave propagation is 
chiefly concerned with the measurement of macroscopic 
elastic moduli, with little interest directed to small scale 
phenomena. After extracting the wave velocity from the 



"coherent" part of the signal, using an appropriate pro- 
cedure, as discussed e.g. in |2^|23, the rest is usually 
discarded as "scattering" or "near field" effects. Most of 
these experiments are done in sands, rather than assem- 
blies of spherical balls (see, however, 25]). 

To summarize: it appears that granular systems can 
be considered as an effective medium for the transmission 
of a short acoustic pulse, if probed on sufficiently large 
scales and pressures, and if one focuses on the initial part 
of the response. This wave front is however followed by 
a noisy tail, which is sensitive to packing details, and 
any quantity or measurement that is dominated by the 
noisy part, such as Cmax.rosp or Cgroup, will not show the 
effective medium behavior. 

Also, probed on smaller scales or possibly at smaller 
pressures, the effective medium description appears to 
become less accurate, even for the initial coherent wave 
front. The two outstanding questions are thus to identify 
under what conditions continuum descriptions hold, and 
if they fail, what other mechanisms come into play. 



III. NUMERICAL MODEL 

Although we performed numerical simulations on 2D 
and 3D granular packings, this paper mainly contains 2D 
results. In our setup (most similar to the experiments of 
Jia et al. plH fT^'l we have a rectangular box containing 
confined spheres under pressure. We send acoustic waves 
through one side of the box and detect force variations 
at the opposite side. 

First we prepare a static configuration of grains. We 
start with a rectangular box filled with a loose granular 
gas. For the 2D simulations the spheres have polydisper- 
sity to avoid crystallization: the diameters are uniformly 
distributed between 0.8 and 1.2 times their average. The 
bottom of the box is a solid wall, we have periodic bound- 
ary conditions on the sides, and the top is a movable 
piston. We apply a fixed force on the (massive) piston, 
introduce a Hertz-Mindlin force law, friction included, 
with some dissipation for the intergrain collisions (see 
appendix . and let the system evolve until all motion 
stops. Our packings are considered to have converged to 
mechanical equilibrium when all grains have acceleration 
less than 10~^° in our reduced units (defined immediately 
below). At this point we have a static granular system 
under external pressure. 

In the rest of the paper we use the following conven- 
tions. Our unit of length is the average grain diameter. 
The unit of mass is set by asserting that the material 
of the grains has unit density. We set the individual 
grain's modified Young modulus E* = 1, which becomes 
the pressure unit (see appendix 0. Since the grains are 
always 3D spheres, we measure pressure even in the 2D 
case as "3D pressure" : force divided by area, where in 2D 
the area is length of box side times grain diameter. The 
speed of sound of the pressure waves inside the grains be- 
comes unity (for zero Poisson ratio). This sets our unit 
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of time, which is about an order of magnitude shorter 
than the time scale of typical granular vibrations. 

Most results are obtained on series of (approximately) 
square 2D samples containing 600 spheres (their centers 
being confined to a plane), prepared with friction coeffi- 
cient /i = 0.5 and Poisson ratio v ~ (see appendix 
The confining stress p that is controlled in the prepara- 
tion procedure, and referred to as the pressure through- 
out this article, is actually the (principal) stress compo- 
nent (T2 (or (T22 in a system of axes for which coordinate 
labeled "2" varies orthogonally to the top and bottom 
solid walls). To check for the influence of p on the re- 
sults, we prepared samples under pressures p = 10~^, 
IQ-^, 10-^ and 10"'' (30 samples for each value). To 
gain statistical accuracy we also prepared an additional 
series of 1000 samples at p = 10"'*. If the particles are 
glass beads this corresponds to 7kPa < p < 7MPa, an 
interval containing the pressure range within which solid 
granular packings are usually probed in static or sound 
propagation experiments. It is worth pointing out that 
the assembling procedure is repeated for each value of 
the pressure. It results in a specific anisotropic equilib- 
rium state of the granular assembly, as characterized in 
section liVAl 

To check the robustness of qualitative results on sound 
propagation, we also studied a few 3D samples, obtained 
similarly to 2D samples by a compaction of a monodis- 
perse granular gas with a piston compressing in the z 
direction to p = a^s = 10"**, and periodic in x and y 
directions. 

In most of the work below, with the exception of sec- 
tion llV F1 we study small amplitude oscillations. For this 
we can use a system linearized around its equilibrium: in 
the static packing we replace the intergrain contacts with 
linear springs, with stiffness obtained from the differen- 
tial stiffnesses (essentially dFt/dt and dF^/dn) of the in- 
dividual Hertz-Mindlin contacts (see appendix 0. The 
equation of motion becomes 

Mii^-Du, (1) 

where for N grains the vector u contains the 3A^ coordi- 
nates and angles (in 2D) of the particle centers, M is a 
diagonal matrix containing grain masses and moments of 
inertia, and D is the dynamical matrix containing infor- 
mation about contact stiffnesses and the network topol- 
ogy. 

Then we solve the eigenproblem of the linear spring 
system and write the oscillation of the grains as the su- 
perposition 

u{t) = a„u(") sin(w„i) . (2) 

n 

The amplitudes a„ are obtained by the projection of the 
initial condition onto the eigenmodes. When the force on 
the top wall is calculated, we calculate the coupling 6„ of 
the given mode with the wall, so the force is 

-Ftop = ^ CLnbn sin(w„i) . (3) 

n 



Typically we will look at the transmission of a short 
pulse through the granular packing. We send in a delta- 
pulse, corresponding to wide band excitation. This cor- 
responds to the following initial condition: at t = the 
grains in contact with the bottom wall have a velocity 
proportional to the stiffness of the contact with the wall. 
This is equivalent to an infinitesimally short square pulse: 
raising the bottom wall for an infinitesimally short time 
and lowering it back. The amplitude of the pulse does not 
matter as all the calculations are linear (except section 
IIV F|) . The quantity we wish to study are the resulting 
force variations on the top wall (the "signal"), defined as 
the force exerted by the vibrating grains on the top wall 
minus the static equilibrium force: 

i^sig(i) =Ftop(i)-^^top(0-) . (4) 

IV. NUMERICAL RESULTS 

In this section we will present the results of extensive 
numerical simulations of sound propagation in granular 
media. We first characterize the geometric state of the 
packing and measure its static properties, including the 
tensor of elastic moduli. Then we report wave propaga- 
tion simulations, showing that in our system, like in the 
experiments of Jia et al. 0, 0| , Fgig is composed of a 
coherent initial peak and an incoherent disordered tail. 
We find that force lines are fairly irrelevant for the pulse 
propagation, and we will show that the decay and spread- 
ing of the initial peak follow packing-detail-independent 
scaling laws. Furthermore we study the pressure depen- 
dence of the transmission velocity of the initial pulse, 
which is compared to the one deduced from static elastic 
properties. We close this section by a short discussion of 
the spectrum and eigenmodes that determine the sound 
propagation, and briefly study the effect of dissipation. 



A. Statics 

We first characterize the system by its geometric and 
static properties, a necessary step as sound propagation 
is sensitive to the internal state of a granular packing. 
Our preparation procedure yields values given in tabled 
for solid fractions $2 (2D definition) and $3 (3D defi- 
nition, in a one-diameter thick layer) and proportion of 
rattlers (particles that transmit no force) /q. For sound 
propagation, the relevant mass density, with our choice 
of units, is equal to the 3D solid fraction $3 of non-rattler 
grains. Values for those parameters are given on account- 
ing for boundary effects: measurements exclude top and 
bottom layers of thickness I, and we checked that results 
became /-independent for I > 3. Tabled also gives the 
bulk values (corrected for wall effects) of the mechan- 
ical coordination number z* (i.e., the average number 
of force-carrying contacts per non-rattler particle), the 
average normal contact force N, normalized by the pres- 
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V 


$2 $3 $3 fo (%) {N/p) Z{2) Z{3) Fa Fi 


10-^ 


0.818 ±0.005 0.561 ±0.005 0.48 ±0.01 3.18 ±0.03 15 ± 2 1.31 1.49 2.81 0.527 0.557 




0.814 ±0.004 0.558 ±0.005 0.49 ±0.01 3.23 ±0.04 14 ± 2 1.27 1.50 2.83 0.524 0.560 


10-5 


0.809 ±0.004 0.554 ±0.005 0.50 ± 0.01 3.31 ± 0.03 11 ±2 1.22 1.49 2.77 0.523 0.558 


10-4 


0.815 ±0.005 0.558 ±0.005 0.52 ± 0.01 3.48 ± 0.03 7.4 ± 1.5 1.06 1.51 2.90 0.524 0.578 



TABLE I: Static data, corrected for boundary effects: solid fractions, coordination number z* , proportion of rattlers fo, average 
normal force, reduced moments (defined in 0) of the force distribution, and fabric parameters F2 (all contacts) and F2 (strong 
contacts) defined in @. Starred quantities are evaluated on discarding rattlers. Stated error bars correspond, here as in all 
subsequent tables, to one standard deviation on each side of the mean. 



sure, and some reduced moments of the distribution of 
normal contact forces, defined as: 



Z{a) 



(5) 



Z values are characteristic of the shape of the force dis- 
tribution (the wider the distribution, the farther from 1 
Z{a) for any a ^ 1). 

The renewal of compression procedure from a granular 
gas at each value of p is responsible for the apparent 
absence of systematic density increases over the range of 
pressure studied here. 

As observed in other studies (see, e.g., 0,113), a di- 
rect compression of a granular gas with friction yields 
rather loose samples, with a low coordination number. 
The minimum value of z* deduced from the condition 
that the number of unknown contact forces should be at 
least equal to the number of equilibrium equations for 
active particles is 3 in 2D, and the values given in ta- 
ble are only barely larger for the smaller values of p. 
This means that our samples have a low degree of force 
indeterminacy at p = 10^^. Complete absence of force 
indeterminacy is obtained with frictionless beads in the 
limit of zero pressure |^ [s^, IH, I3 ■ With frictional con- 
tacts, it appears that some force indeterminacy persists 
even in the limit of zero pressure [sO, IH, [s^ . 

The fraction of rattlers is very large (up to 15%), when 
z* is close to 3, while the number of rattlers fastly de- 
crease as z* increases with p. We also computed the state 
of stress in the samples: due to the assembling proce- 
dure, the principal directions are horizontal (coordinate 
X, principal value ai) and vertical (coordinate y, princi- 
pal value 172), with a larger vertical stress, Gxjai < 1. 
This ratio (table ITl|l stays constant, within statistical er- 
ror bars, as the controlled confining stress a2 = p is var- 
ied, except at the highest pressure p = IQ-"*. (Stresses 
are readily evaluated with the usual formula 



y Z^^v' 



P 



where is the a component of the force exerted on par- 
ticle i by particle j, and the vector rij, pointing from the 
center of particle i to the center of j, should account for 
periodic boundary conditions and involve "nearest im- 
age" neighbors). 



The force network of a typical granular configuration 
of 600 grains prepared at pressure p = lO^* and friction 
coefficient = 0.5 is shown on Fig. QJa). We show the 
stiffness network on Fig.^b) and the histogram of forces 
and stiffnesses on Fig. Q^c), to which we will come back 
shortly. It is interesting to note that the shape of the 
force distribution changes very little in the investigated 
pressure range: values of reduced moments Z, listed in 
table U hardly change with pressure. In ref. [S^, a "par- 
ticipation number" 11 was defined, as an indicator of the 
width of the force distribution, or the "degree of localiza- 
tion" of stresses on force chains. In fact, one would have 
n = \/Z{2) if we had defined Z with the magnitude of 
the total contact force instead of its sole normal compo- 
nent. Makse et al. js^l found that 11 increases linearly 
with Inp, until it saturates ai p > 10-'"'. Our contra- 
dictory observation of a nearly constant Z{2) is likely 
due to our compressing the system anew from a granular 
gas at each p, instead of quasistatically increasing p in a 
previously assembled solid sample. 

The anisotropy of the contact network (which carries 
anisotropic stresses) is apparent in the histogram of con- 
tact angles and the force histograms. Figure |21 shows the 
histogram of contact directions. We also plot the con- 
tribution from strong (contact force larger than the av- 
erage) and weak (force smaller than average) contacts. 
Due to the assembly procedure there is an anistropic 
stress field, resulting in a bias of the strong contacts to- 
wards the vertical direction. Let us recall that this is also 
the principal direction corresponding to the largest eigen- 
value of the stress tensor. As a direct consequence of the 
force anisotropy, the small forces are biased in the oppo- 
site direction, although this effect is less pronounced |36| . 
One way to quantitatively assess the importance of this 
anisotropy is to compute the fabric parameter 



(6) 



where Uy is the vertical coordinate of the unit normal 
vector and the average runs over the contacts. The de- 
parture of F2 from its isotropic value 1/2 measures the 
anisotropy. Table Ogives F2 for the different pressure lev- 
els, along with parameter F2 , obtained on counting only 
the contacts that carry larger than average forces. Once 
again, the level of anisotropy does not depend on pres- 
sure, except for a slight difference at p = 10"^, in which 
case it is a little larger (consistently with the larger stress 
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FIG. 1: (a) Snapshot of a force network at pressure p = 10"*. Line widths are proportional to the force between the grains, 
grain centers are plotted with gray dots. Only the normal forces are plotted, the tangential (frictional) are not. (b) The stiffness 
network of the same configuration. Line widths are proportional to the stiffness dFn/dn of the contact (normal part). While 
the force network shows considerable spatial fluctuations, the stiffness network is much more homogeneous, (c) Histogram of 
the normal contact forces and contact stiffnesses of fOOO configurations. The area under the two curves are the same. The plot 
shows that the stiffness is more narrowly distributed, a feature discussed in more detail in Sec. ITVni 



all contacts strong contacts weak contacts 

FIG. 2: Histogram of contact directions at pressure p = 10~*. The diagrams show the average of 1000 independent con- 
figurations, 600 grains each, and only bulk contacts are counted (at least 3 particle diameters away from walls). The strong 
contacts (contact force larger than the average) show significant anisotropy: vertical directions are favored. The weak forces 
are much more isotropic, although slightly more of them are horizontal than vertical (as described in [s^). Recall that the 
piston compressed the grains in vertical direction. 



anisotropy). 

If, on long length scales, the material can be consid- 
ered as an ordinary homogeneous 2D material with two 
orthogonal symmetry axes, the granular packing has 4 in- 
dependent macroscopic elastic moduli, which relate the 
3 independent coordinates of the symmetric stress tensor 
(Ty to the 3 coordinates of the symmetric strain tensor 
Ci i as 



u as 



o-ii 








0'12 





Cll Cl2 
Ci2 C22 








2C33 





£11 




£22 




_£12_ 



(7) 



In Eqn. [7| indices 1 and 2 correspond to the horizontal 
(periodic) and vertical (along which the normal stress is 
controlled) directions on the figures. Counting positively 
shrinking deformations and compressive stresses, strain 
components should be related to the displacement field 



duj 



9uj 
dx. 



Principal stress ratios a\l 02 and apparent values of the 4 
moduli introduced in Eqn.[7| obtained on imposing global 
strains on the rectangular cell containing the samples, 
are given in table |n] To check for possible length scale 
effects on static elasticity (in view of the small system 
size), samples were submitted to inhomogeneous force 
fields, or to various local conditions on displacements. 
The result of these computer experiments, described in 
Appendix^ is that constants C22 and C33 are already 
quite well defined on the (modest) scale of the 600 sphere 
samples (typically 24 x 25). The assumption that non- 
uniform stress and strain fields, which vary on the scale 
of a fraction of the sample size, are related by iQ, with 
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Pressure 


0"l/o"2 Cll 


Cl2 


C22 




C33 


10-^ 


0.79 ±0.06 (9.9 ± 0.9) X 10-"* 


(8.1 ±0.4) X 10-4 


(13.9 ±0.8) 


X 10-4 


(1.2 ±0.4) X 10-4 




0.79 ±0.06 (2.3 ± 0.2) X 10"^ 


(1.7 ±0.08) X 10-^ 


(3.2 ±0.2) 


X 10-^ 


(0.35 ±0.07) X 10-^ 


10"^ 


0.78 ±0.06 (5.4 ± 0.5) X 10-^ 


(3.5 ±0.2) X 10-=* 


(7.4 ±0.3) 


X 10-^ 


(1.1 ±0.2) X 10-^ 


10-4 


0.71 ± 0.04 (1.35 ± 0.09) x 10-^ 


(6.8 ±0.4) X 10-^ 


(1.88 ±0.08 


) X 10-2 


(4.0 ±0.5) X 10-^ 



TABLE 11: Stress ratio and elastic moduli, as defined in Q, for the 4 investigated pressures with ^ — 0.5 and p — 0. 



the elastic constants of table ^ predicts results which 
approximately agree with numerical tests on our discrete 
packings, in spite of their moderate size. This agreement 
is better for the longitudinal constant C22 than for the 
shear modulus C33, and improves on increasing p. 

The conclusion that macroscopic elasticity applies even 
at moderate length scales is in agreement with the results 
by Goldhirsch and Goldenberg on homogeneously forced 
disordered packings [s^- However, when probing the re- 
sponse to localized forces (which perturb the system inho- 
mogenously even in the elastic limit) those authors identi- 
fied a larger length scale of about 100 diameters in order 
to recover macroscopic elasticity [SS]; these differences 
might also be due to the fact that their study concerned 
frictionless quasi-ordered systems. We should keep in 
mind also that Goldhirsch and Goldenberg looked at the 
full spatial dependence of the elastic response, while we 
extracted only global elastic quantities. One may also 
note (see, e.g., |39|) that constitutive laws are obtained 
with numerical simulations of disordered granular sam- 
ples in the quasistatic regime with relatively small finite 
size effects when the number of particles is above 1000, 
and that the level of uncertainty and fluctuations is fur- 
ther reduced on investigating the response to small per- 
turbations on a fixed contact network. Tanguy et al. ^3 
studied the finite-size effects on the elastic properties of 
2D Lennard- Jones systems at zero temperature. While 
large sample sizes (N ~ 10000 particles) were necessary 
for the low-frequency eigenmodes to resemble the macro- 
scopic predictions (an issue we shall return to in the se- 
quel to this paper) they did observe apparent elastic con- 
stants, as measured on globally deforming the sample, 
to converge quickly (for N ~ 100) to their macroscopic 
limit. The observation of macroscopic elastic behavior 
(with some limited accuracy) in an assembly of 600 par- 
ticles is not really surprising in this context. 

Assuming macroscopic elasticity to hold in our sam- 
ples, elastic constants C22 and C33 determine velocities 
of longitudinal (q) and transverse (ct) waves propagat- 
ing in direction 2 (normal to top and bottom walls), as 



ce 



andct = J— . 



(8) 



It should be pointed that the different elastic constants 
do not exhibit the same scaling with the pressure. Most 
notably, the ratio of shear modulus C33 to longitudinal 
modulus C22 (or the ratio of the corresponding wave ve- 





FIG. 3: Snapshot of the oscillations. The length of the arrows 
show the (magnified) displacement of the grains from their 
equilibrium position. One can see the localized large ampli- 
tude oscillations of the grains near the bottom wall (source), 
and a smaller amplitude homogeneous wave traveling towards 
the top wall (detector). At the time of the snapshot, t — 80, 
the wave almost reached the top wall. 



locities, according to steadily increases with p. We 
shall return to this issue, and compare different predic- 
tions for sound velocities and their pressure dependence, 
in section HV Dl and in the discussion. 



B. Wave propagation: qualitative observations 

We now turn our attention to the pulse propagation. 
Before studying -Fsig, we will discuss here an example 
of the spatial structure of the propagating pulse. Our 
first observation is that acoustical waves do not correlate 
in any obvious way with the existence of force chain- 
like configurations. A snapshot of the grain oscillations 
shortly after the system is "kicked" is shown on Fig. O 

This figure clearly shows that the naive idea that the 
acoustic waves would follow the strongest granular force 
chains is false. Instead, one can see the propagation of a 
rough wave front. One reason we can immediately point 
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FIG. 4: (Color online) The signal -Fsig, the extra force exerted 
by moving grains on the top wall, (a) Four independent 2D 
frictional configurations are shown, the signal corresponding 
to the packing of Fig. Qis plotted in red. The arrow indicates 
the time of the snapshot of Fig. |3 (b) The ensemble average 
and root-mean-square deviation of 1000 independent configu- 
rations. The first cycle of the oscillations is almost the same 
on all configurations (we call this the coherent part of the 
signal), while the following part is very much configuration 
dependent. The ensemble average only contains the coherent 
part plus some weak broadened sign of multiple refiections on 
the top and bottom walls, (c) 2D frictionless and (d) 3D 
frictionless systems exhibit similar behavior. Time (as well 
as length in subsequent figures) is denoted in dimensionless 
units, see Sec. Illll for details. 



to is that even though the forces of the intergrain con- 
tacts exhibit a strong spatial fluctuation, the stiffness is 
much more homogeneous, see Fig.^b). This can be un- 
derstood simply as follows. Consider a force law which 
in scaled units reads — n^, where n is the normal 
deformation. The stiffness s is then simply given by 
dFn/dn — f3n^~'^. Clearly, for /3 = 1 (corresponding 
to the 2D Hertzian force law), all the stiffnesses values 
are the same. For the Hertz-Mindlin law, /? = 3/2, and 
we find that the stiffnesses are proportional to the cubic 
root of the contact forces, leading to the rather homo- 
geneous stiffness network shown in Fig. ^b). So if we 
compare two links with forces differing by a factor 8, the 
corresponding stiffnesses only differ by a factor 2, and 
the sound speed — proportional to the square root of the 
stiffness — differ only by a factor of \/2- Even though the 
contact forces follow a wide distribution, the stiffness- 
distribution is strongly peaked [see Fig.^c)]. Although 
this is a rather trivial observation for Hertzian contacts, 
we are not aware of its being explicitly mentioned in the 
literature. 

An additional reason for the weak effect of force chains 
on the sound propagation may be that the disorder of the 
grains is significant: on a force chain with weak side links 
the oscillation quickly spreads into its neighborhood, re- 
sulting in a more homogeneous base of the oscillations. 
Anyway, the conclusion we can draw here is that the 
force chains are not relevant for the evolution of the ini- 
tial wavefront. 



C. The coherent wavefront 

Let us now study the experimentally accessible signal 
-Fsig- The time dependence of this signal is shown on 
Fig. E^a). Clearly Fsig can be thought to be composed 
of an initial peak followed by a long incoherent tail. One 
can see that for configurations that are similar in overall 
geometry but statistically independent, the initial first 
cycle of the signal is very similar, but the following part is 
strongly configuration dependent. The time dependence 
of the signal is very reminiscent of the traces measured by 
Jia et al. |T5l | in their ultrasound experiments. Following 
the nomenclature introduced by these authors, we call 
the first part of the signal the coherent part. In the 
ensemble average only the coherent part of the signal 
shows up (plus its later weaker echoes) [see Fig. Elb)]. 
The random part of the signal contributes to the root- 
mean-square deviation. We also found that qualitatively, 
Fsig is very similar for 2D frictional, 2D frictionless and 
3D frictionless systems, see Fig. 01[a,c,d). 

We will now focus on the initial peak of Fsig, and de- 
termine for this coherent wavefront its propagation ve- 
locity, and the time evolution of its shape. We have only 
measured the time dependence of the signal at a fixed 
distance, and a qualitative picture of the evolution of the 
coherent wave can be extracted from a sequence of mea- 
surements at varying source-detector distance. This is 
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FIG. 5: The coherent part of the signal in containers of 
varying height (source-detector distance). For taller contain- 
ers the signal arrives later with decreased amplitude and in- 
creased width. These are quantitatively analyzed in the next 
few figures. 



shown on Fig. during the propagation of the signal (as 
it arrives later at longer distances), the coherent part's 
amplitude decreases, and its width increases. 

Now we look at Fgig quantitatively. As shown in the 
inset of Fig.|Bl we characterize the coherent peak by three 
points: its peak location, its first 10% of peak value and 
its first zero crossing. In Fig. we show, for various 
source-detector distances, the times at which these three 
characteristics can be observed at the detector. In rea- 
sonable approximation, the time of flight depends lin- 
early on the source-detector distance, although the up- 
ward curving of the data suggest that for small systems 
the propagation velocity appears larger than for large 
systems. We define the time of flight velocity, Ctof, by 
measuring the difference of the arrival times (at 10% of 
the peak's level) for source-detector distances of 8 and 25. 
The velocity thus defined can be measured in reasonable 
small systems (containing 200 and 600 particles respec- 
tively), while on the other hand being quite close (within 
10%) from the large scale velocity. Based on this defini- 
tion of time of flight, we have a sound speed, Ctof = 0.26 
in our units, at pressure p — 10~*. 

In Fig.[7|we plot the scaling of the amplitude and the 
width of the coherent part of the signal. The amplitude 
is well approximated with a power law, A ~ for the 
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FIG. 6: The arrival time of the coherent part of the signal as 
a function of the source-detector distance. The inset shows 
the definition of the symbols: leading edge at 10% of the 
first peak height (A), the first peak (0)i ''■nd the first zero 
crossing of the signal (□). All three characteristic points of 
the signal have a linear time-distance relation. The slope of 
the time-distance plot of the leading edge defines a time-of- 
flight velocity: Ctof = 0.25. 



2D simulations 7 « 1.5. The width of the coherent part 
of the signal increases with distance also as a power law: 
~ L". For the 2D simulations the increase is close to 
linear, awl. 

We are not aware of any prediction or previous analy- 
sis of these exponents 7 and a for polydisperse random 
packings. In order to put these results into perspective, 
it is important to keep in mind that i^sig is not the ampli- 
tude of the wave motion in the medium, but the resulting 
force on the boundary at the other edge. Since the force 
is proportional to the local stretching, i.e., the deriva- 
tive of the amplitude of the wave, 7 is not the exponent 
with which the wave amplitude itself decays [see also the 
discussion at Eq. lfTB)l ]. 

We have compared this behavior with the behavior of 
propagating pulses in a one-dimensional chain of balls. 
Even in this simple system, dispersion effects (wavenum- 
ber dependence of the frequency of the waves) give rise 
to nontrivial exponents — as we shall discuss in more de- 
tail in Sec. |Vl both the exponent and the shape of the 
pulse can be determined analytically. We collect the ex- 
ponents 7 and a in Table ITTll An important lesson from 
the ID analysis is that the decay exponent 7 is not uni- 
versal, as it depends on the precise shape of the initial 
pulse: 7 = 2/3 for our usual initial condition (equilib- 
rium position but nonzero velocity next to the wall at 
t — 0), and 7 = 1 if the initial condition is zero velocity 
but nonzero displacement at the wall (not plotted). If we 
allow polydispersity in the ID chain, the scaling appears 
to have a larger exponent depending on the magnitude 
of the polydispersity, although from the data shown in 
Fig. [71 we cannot draw a definite conclusion. 

In conclusion, the main qualitative differences between 
the ID results and those for the coherent pulse in the 
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TABLE III: The scaling exponents 7 and a for different gran- 
ular systems. 



FIG. 7: The scaling of the amplitude and width of the co- 
herent part of the signal with the source-detector distance 
L. Upper panel: the amplitude follows roughly A ^ 
(inset). In the main panel we plot the effective value of the 
exponent: 7cff = dlog A/d\og L. The symbols are 2D disor- 
dered (full circles), ID chain of identical balls (open circles) 
and ID chain of polydisperse balls (other symbols, with vary- 
ing polydispersity) . Lower panel: the width of the coher- 
ent part of the signal increases with distance: width --^ L°' 
(inset). In the main panel here also the effective exponent 
Ooff = dlog W/ dlog L is plotted. 



disordered 2D packings is that (i) in the ID chain the first 
pulse broadens as t^^^ whereas the pulse in the disordered 
2D medium broadens linearly; (ii) the amplitude of the 
pulse decays much faster in the disordered medium than 
in the ID chains (in other words, 7 is larger). 



D. Speed of sound, elastic moduli, and pressure 
dependence 

In this section we turn our attention to the sound 
speed, and in particular study its variation as a func- 
tion of the confining pressure p. The main quantity is 
the time-of-flight velocity obtained from the propagation 
of the coherent pulse Ctof (see section ITV C|) . It should 
be compared to the values of transversal (ct) and longi- 
tudinal (q) wave speeds that are deduced (Eqn.lHJ from 
the apparent elastic moduli of tabled We also compare 
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our results to experimental data for sound propagation; 
since some experiments have been performed on regular 
packings, we also have studied these analytically and nu- 
merically (see section An overview of these various 
propagation velocities as function of pressure is shown in 
Fig. IHl Let us first discuss the scaling of ctof, ce and ct 
(as defined in 0) with p. Recall that for a fixed contact 
network with Hertzian forces which stay proportional to 
the pressure p, the sound velocity scales as p^^^. We find 
here that Ctof follows this scaling quite accurately, while 
Q appears to be growing slightly faster as p^'^^. Surpris- 
ingly, data for the velocity ct of transverse waves abide 



.0.23. 



we do not know the 



by a different scaling, ct ~ p' 
reason for this behavior. 

Since the coherent wave is essentially longitudinal in 
nature, one should compare ci and ctof- Even though 
both quantities scale rather similarly, ctof is roughly 40% 
larger than q. As discussed in section ITV CI our defini- 
tion of Ctof is based on measurements in relatively small 
systems, and from a few simulations in larger systems we 
found that this may overestimate ctof by some 10%. In 
addition, if we do not measure the first arrival of the sig- 
nal, but instead measure the first peak location, or the 
first zero crossing, ctof would go down substantially. 

Furthermore, it seems that the pulse propagation with 
our method of excitation does not probe the material on 
the longest scale. On shorter scales, the material appears 
somewhat stiffen As discussed in appendix^ numerical 
measurements of the elastic modulus C22 can be per- 
formed on various length scales, the shorter ones, in the 
case when displacements are locally controlled, leading 
to larger apparent values of C22- In addition, we shall 
show in section llV El that there is a strong contribution 
from non-plane- wave modes which cannot be expected to 
be described by continuum elasticity. 

It might therefore be concluded that a simple long- 
wavelength description gives a good first approximation 
of the propagation velocity of the coherent wavefront, 
but that modes that are not accurately described by a 
long-wavelength approximation contribute substantially 
to the wave propagation for the system sizes and exci- 
tation method employed here. For a triangular lattice 
of monodisperse balls, we compare the analytical expres- 
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FIG. 8: (Color online) The pressure dependence of the various sound speeds. The main data are the results for Ctof obtained 
from our simulations, which show perfect p^^® scaling. Velocities C£ and Ct deduced from elastic moduli, as in JSJ, are smaller; 
one would have expected that Ctof ~ Q. Ct is much smaller and has been multiplied by 2 to fit within the scale of the plot. The 
theoretical curves for triangular lattice are Eqs. I18II (frictionless) and 1191 (frictional). For the latter there is a slight variation 
depending on the Poisson ratio of the grain's material: the gray band corresponds to the range < v < 0.5, Simulations 
for the frictionless triangular lattice (+) show excellent agreement with Eq. 1181 . For the frictional case (x) the simulation 
shows significant finite size scaling: it should approach the top side of the gray band (we used i/ — and size L = 24, but for 
p = 5 X 10~^ a larger system L = 160 is also plotted). The simulation for 2D disordered frictional case (O) shows results very 
similar to the triangular lattice. For comparison we also show three experimental datasets [recall the pressure- and velocity 
scales: E* = E/{1 — u^) and c* = E* / p ]: triangular lattice of steel spheres (green □, from [l^). triangular lattice of nylon 
balls (red Qi S'lso from p^: see text for explanation of the arrows), and disordered 3D glass spheres (blue A from 15]). For 
reference, lines with slope 1/6 (for the p^^^ law) and 1/4 (sometimes quoted as effective exponent for low pressures) are shown. 



sions for the sound speed Eq. (|18|l and H19|l . which are 
derived in section^for infinitely large lattices, with sim- 
ulations on finite lattices. Both the frictionless and fric- 
tional cases are in excellent agreement, even though the 
frictional one shows appreciable finite size corrections. 
The simulation for 2D disordered frictional case shows 
results very similar to the triangular lattice — including 
the p^^^ scaling expected naively from the Hertzian force 
law — for the range of pressures considered. 

This quantitative agreement is somewhat surprising in 
view of the large difference in coordination numbers (6 
versus barely larger than 3). This might partly be due 
to the small wavelength effects, which affect the results 
in disordered systems, while one easily observes the long 
wavelength result Ctof = with a perfect regular lattice. 

The only 2D data we are aware of are for spheres on a 
triangular lattice. These systems are inevitably slightly 
polydisperse, which prevents the closing of all contacts 
between nearest neighbors on the lattice 3, 17, and 
in the limit of low pressure, the coordination number 
should not exceed 4. However, once the reduced pressure 
is high enough for the elastic deflection of contacts to 
compensate for the open gaps, the behavior of the perfect 
lattice is retrieved. This effect can be evaluated with the 



reduced pressure defined in [41| as 

where ad is the width of the diameter distribution. Ef- 
fects of polydispersity disappear as P* grows beyond 1. 
For steel spheres the data of |OS (for which a ^ 10""*) 
fall close to our calculated values for the triangular lattice 
with friction when p > 10~^, as expected. Even though 
there is a discrepancy in the velocity of the order of 10- 
20%, this agreement is remarkable, since Ctof has been 
calculated without any adjustable parameters. Possible 
finite size effects might explain why these data lie below 
the theoretical frictional curve for the perfect lattice. 

The triangular lattice of nylon balls |3 shows signif- 
icantly larger rescaled velocity than expected. Possibly, 
this discrepancy is simply a reflection of the uncertainty 
in the effective elastic constant at the frequency range 
of the experiments: nylon is a viscoelastic material for 
which the Young modulus increases strongly with fre- 
quency. We do not know the values of the elastic con- 
stants at the experiment's frequencies, but nevertheless 
if we use a Young modulus twice as large as its zero fre- 
quency value (for the plot the zero frequency modulus 
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FIG. 9: (a) Eigenfrequencies of the linear system for the 
packing shown on Fig. The squared eigenfrequencies are 
plotted against the number of the mode n. Modes n — 
0, . . . 140 have eigenvalue zero, as a consequence of "rattler" 
grains which are not connected to the network. The inset 
shows a magnification of the plot around the first few nonzero 
eigenvalues, (b) The contribution of the eigenmodes to the 
transmission signal anb„ (see Eq. O. On both panels the 
eigenmodes plotted on Fig. llOl are marked by circles. 



was used), then the curve would shift as indicated by the 
arrows. 

Finally, disordered 3D glass spheres 15j display smaller 
velocities than any 2D case. One possible explanation 
is that 2D experiments on planar sphere assemblies can 
be viewed, if we imagine stacking such layers on top of 
one another, as probing the stiffness or wave propaga- 
tion along dense, well coordinated planes in a 3D mate- 
rial with extreme anisotropy. This renders plausible the 
observation of unusually high sound velocities, in com- 
parison with ordinary 3D packings. 



E. Eigenfrequencies and eigenmodes 

In a rectangular sample of homogeneous elastic mate- 
rial with boundary conditions similar to those employed 
here, the eigenmodes are plane waves with wavevector 
k = (fci,fc2) — {nij^,n2^). If the tensor of elastic 
moduli has the form given in Q, then it is straightfor- 
ward to show that the associated frequencies u;+ , w_ are 






FIG. 10: A few selected eigenmodes of the linear system, 
(a) n = 141 and (b) n — 142 are the first two nonzero eigen- 
modes. They correspond to the lowest excitation modes of a 
continuum body, though slightly distorted by the disordered 
contact network, (c) n — 197, (d) n = 674, and (e) n = 974 
are some of the modes that contribute significantly to the 
transmission of the signal, (f) n = 1707 is a high frequency 
localized mode. The modes shown here are marked on the 
eigenvalue plot, Fig.^ 



given by = A±/$^, 
acoustic tensor 



A± being the eigenvalues of the 



A(fci,A:2) 



Ciifc? + C33fc| (Ci2 + C33)hk2 
{Ci2 + C33)fclfe C33fc? + C22k^ 



which implies that uj on k in the long wavelength limit. 

We show the spectrum of eigenmodes for a granular 
packing of 600 grains at pressure 10^^ in Fig.El^a), and a 
few selected eigenmodes in Fig.^1 There are a number of 
zero eigenvalues because of "rattler" grains not connected 
to the force network. The lowest nonzero modes corre- 
spond to (slightly distorted) solid body modes, which are 
similar to those expected from continuum theory. Re- 
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FIG. 11: The transmitted signal with damping. The level of 
damping is expressed as a fraction of the critical damping on 
each contact. The damping affects the coherent part of the 
signal much less than the random part. 



markably, in the absence of friction it is much harder to 
identify eigenmodes corresponding to continuum media 
modes, even for low frequencies; and the low frequency 
modes are more abundant (not shown on the figures). 
Nevertheless, the transmission signal looks rather similar 
to the frictional case (Fig.0J. 

There are a large number of localized eigenmodes 
[Fig. liur f)]. which do not contribute substantially to the 
signal transmission; clearly the modes that dominate the 
transmission are global modes [Fig. I^b)]: they contain 
oscillating grains at both at the source and the detector 
wall. But with the exception of only a few modes (with 
mode number 141-147 roughly), their appearance is quite 
different from simple plane waves [Fig. Iior c-e)]. This 
indicates that, at least for the system sizes, pressures 
and excitation method employed here, the transmission 
of sound cannot be captured completely by considering 
the material as a simple bulk elastic material. In fact, in 
the light of these findings it is remarkable how close the 
continuum prediction C£ comes to Ctof. We will present a 
more extensive study on these eigenmodes elsewhere |42j . 



F. The effects of damping 

Finally, we show here how damping affects the wave 
propagation. We added viscous dissipation to the Hertz 
contacts in the way described in appendix ^ The re- 
sulting system cannot be easily described by a linear 
model, and we obtained the wave propagation signal by 
molecular dynamics simulations of the grain oscillations. 
On Fig. ^] we show the transmission signal for a single 
configuration with various levels of damping. For large 
damping the coherent part of the signal is only slightly 
altered, while the random part is strongly suppressed. 
This is in qualitative agreement with the experiments of 
Ref. [T^ . where damping was induced by adding a small 
amount of water to the glassbead packing. 



V. ANALYTIC RESULTS 
A. ID chain 

The problem of the propagation of a pulse in a ID 
granular chain has been considered by many authors 
E El El El El El El but the majority of the work 
is concentrated on analyzing an initially uncompressed 
chain. In this case the nonlinear force law plays an im- 
portant role, as well as the fact that there are no restor- 
ing forces between the balls which initially just touch. 
For comparison with sound propagation in granular me- 
dia as a function of pressure, the relevant approach is 
to first linearize the equations of motion starting from a 
compressed chain, and then study the propagation of the 
pulse as governed by these linearized equations. 

The simplest system resembling the problem of sound 
propagation in granular media (under pressure) is a ID 
chain of identical elastic balls, confined and compressed 
between two walls. At t = we disturb the first ball (see 
below for details), this disturbance travels with sound 
speed c in the chain, and arrives at the other wall at 
time to = N£/c, where £ is the diameter of the balls. For 
this system we can calculate the scaling exponents and 
the waveform analytically. 

In the Appendix we calculate the attenuation exponent 
of this wave. For initial condition (A), where the first ball 
has nonzero velocity but zero displacement at t = 0, the 
force with which the last ball presses the wall at time to 
scales with N as 

F^{to) - 7V-'/^ . (9) 

Initial condition (B), where all balls start with zero veloc- 
ity but the first has a finite displacement, gives a different 
answer: 

F''ito)^N-K (10) 

These are the attenuation exponents for a uniform ID 
chain. 
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To derive the waveform analytically in the large sys- 
tem and long time limit, we consider the long wavelength 
expansion of the dispersion relation ljC2() : 

a;„ « cK - -^kl ■ (11) 

A propagating wave solution u{x,t) — Aexp{ikx — cut), 
where for long wavelengths x can be considered continu- 
ous, has to satisfy the following partial differential equa- 
tion to match dispersion relation Hll|) : 

du du cf' d^u 

"9? """a^ ^ 24 ^ ' 

Changing variables to the co-moving frame, ^ — x — ct^ 
the du/dx term drops out. Looking at similarity solu- 
tions of the form 

u{^,t) ^ t-m [w ^ , (13) 
we obtain a = 1/3 and 

= -gU{w) - |c/'(^) + ^C/"'(u;) ■ (W) 

This leads to different classes of solutions for different at- 
tenuation exponent g. First we consider the case 17 = 0, 
which leads to Airy functions: Uq{'w) = Ai(2w/(c^^)-^/^). 
We can generate further solutions by differentiating 
Eq. p4|l . For example by differentiating once and twice 
we find that U = Uq solves Eq. (tTIjl for g = 1/3, while 
Uq solves it for g — 2/3. For u{x,t) this gives solu- 
tions, e.g., for g = 1/3 as ^-^/^ Ai[2(x - ct)/iceHy^^], 
which as we show soon is the selected solution for initial 
condition (A). At this point we need more information 
to see which solution is selected: Eqs. © and ((117)) tell 
the exponent of the scaling of the force at the wall with 
N. Note that u{x,t), is the propagating solution in a 
semi-infinite medium; the solution for a reflecting wall 
boundary condition, u(x = (A'^-l- l)^,t) = 0, is composed 
of two counter-propagating waves. Now the force at the 
wall is proportional to the displacement of the ball (at 
X = Ne) next to the reflecting wall [at x = (A^ -I- 1)^], 
hence it can be written as 

F{t) = K{u{x, t) - u{2{N + l)e ~x,t)} , (15) 

which with x — N£ and t — to + t gives for the 5 = 1/3 
solution 

^<'« -^-'-^ (7) w^)- 

— («-t)""-'(w^^) 

Note that because of the extra differentiation the decay 
exponent of the force on N does not equal g — 1/3 but 
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FIG. 12: (Color online) Comparison of signal shapes for 
initial condition (A). The theoretical prediction of Eq. H16|l 
(thick black line) is compared to the ID chain of identical 
balls, numerical sum in Eqs. <C4HC6t (color lines). The sim- 
ulation of a perfect triangular lattice (red full circles) is very 
close to the ID chain of the same size. Note that Eq. (1161 has 
an undetermined multiplicative factor as it is obtained as a 
solution of a linear equation. 

instead becomes 7 = 5 + 1/3 = 2/3. This scaling expo- 
nent is the same as for the initial condition (A) in Eq. (|2J), 
showing that indeed the 5 = 1/3 solution is selected here. 

To match the initial condition (B), we need to use the 
5 = 2/3 solution: 

•^--(»-?)"'-'"(l]^^l^)' 

The solution of the ID chain, obtained by numerically 
evaluating the sum (|C4|) converges for A^ — > 00 to the 
analytical solution, see Fig.ll2lwhere the initial condition 
(A) is plotted. 

At this point we can understand the connection be- 
tween the two initial conditions. Initial condition (B) is 
related to (A) by a time differentiation. Since the equa- 
tions are linear, the solutions are similarly related to each 
other. The above solutions have the structure that differ- 
entiating one of them and dropping subdominant terms 
gives another of the solutions, with exponent g which has 
increased by 1/3. 

If we allow for disorder in the ID chain, the results 
change slightly. We introduce disorder by varying the 
radii of the ball, and solve this ID system numerically, 
as described in Sec. IIIII for 2D and 3D packings. On 
Fig. 13 we show the exponents of the polydisperse ID 
chain. Both the amplitude exponent 7 and the width 
exponent a appear to be larger than in the case of iden- 
tical balls, but the results are not clear enough to extract 
a value for the exponents. 

The above analysis shows that in the long time limit, 
the propagation of an initially localized pulse is governed 
by an Airy equation — as Fig. E| shows, the first pulse 
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and the first oscillations behind it converge to an Airy 
function type behavior when viewed in a frame co-moving 
with the initial pulse. Note that the kinetic energy in the 
leading pulse decays rapidly as t^^^^^/'^: 



kin 



init. pulse 



m 
"2" 



init. pulse 



7273" 



t2/3 



^1/3 ^ ^-2g-l/3 



(17) 



because the number of terms in the sum which contribute 
to the first peak is proportional to the width of the pulse, 
which scales as t^/'^. Hence for the pulse shown in Fig. 1121 
the kinetic energy in the first pulse decays as since 
g = 1/3. This illustrates that as time progresses, more 
and more of the energy is stored in the region behind 
the first pulse. The oscillations in this region are rel- 
atively incoherent, with a frequency comparable to the 
maximum frequency of the dispersion relation. As the 
size of this region increases linearly with time, the typi- 
cal amplitude of these oscillations decays as t"^/^. One 
can also obtain a t~^^'^ type decay directly from a steep- 
est descent analysis near the maximum of the dispersion 
relation of the linearized equations of motion. 

The fact that the first pulse in ID chains is described 
by an Airy function has been noted before 0, 23l ■ Most 
of these studies are for initially uncompressed chains, 
however. In this case, all the energy remains confined 
in the first pulse, due to the absence of restoring forces. 
As a result, the exponent of the time dependence of the 
amplitude is different, and consistent with energy con- 
servation in the leading pulse. 



B. Triangular lattice 

In view of the experiments of Gilles and Coste fljl , it 
is illuminating to also apply these results to the simplest 
2D system: a triangular lattice of balls, with rectangu- 
lar boundaries. The initial condition is given on balls 
touching one wall of the rectangle, and we assume that 
a lattice direction of the triangular lattice is parallel to 
this wall. The longitudinal sound speed in a perfect tri- 
angular lattice (no polydispersity) of Hertzian balls can 
be easily calculated, see, e.g., ^j. For the frictionless and 
frictional case respectively it is given by 



□fric 



„fric 



23/27ri/2 \E*) ' 

23/2^1/2 3 



1/6 



(18) 
(19) 



where the parameter rj is the ratio of the tangential and 
radial stiffnesses of a Hertz-Mindlin contact, see Eq. ljA4|l . 

One way to calculate this is to map to the ID chain 
of identical balls. If in the triangular lattice the longitu- 
dinal motion is perpendicular to rows, then a row of M 



balls moves together, corresponding to a single ball in ID. 
Thus of the N rows each has mass ttLcs — Mpn/Q, they 
are separated by distance ^ = \/3/2 (recall our length 
unit was the ball diameter), and connected by an effec- 
tive spring i^eff — 3^/^2~27Vfj)i/3_ \^\yq frictional case 
K^s has an additional factor of (1 + 7^/3). 

This way we also predict the shape of the signal for 
the triangular lattice, see Fig.^l The cause of the slight 
deviation from the ID chain result is a consequence of the 
fact that in the triangular lattice the springs connecting 
to the walls are different: -fCwaiiZ-K'buik = 3^/^/4 w 0.62. 

If the radii of the balls are polydisperse, then at pres- 
sures low enough (that the length scale of the elastic de- 
formations become comparable to the polydispersity) the 
stress field fluctuates spatially. The effect of this on the 
sound speed has been calculated by Velicky and Caroli 
^ in a mean-field approximation. 



VI. DISCUSSION 

We presented numerical simulations of pulse propaga- 
tion in ID, 2D and 3D granular systems. This response 
can be decomposed into an initial coherent part, which is 
independent on the details of the packing, and a subse- 
quent random part, which is strongly realization depen- 
dent. We have focused on the properties of the initial 
coherent front. Our first observation is that the response 
to a pulse propagates linearly in time, defining a time-of- 
flight velocity, and does not follow force chains. 

The fact that the packings in our numerical simulations 
have roughly the same number of grains per container 
side (although in 2D) as the systems which have been 
studied experimentally by Jia et al. p^ . and that our 
temporal signals are very comparable to the experimen- 
tal ones, makes us confident that our simulation results 
can be fruitfully compared to experiments like these. In- 
deed we find that the 3D experimental and 2D numerical 
results for the time-of-fiight velocity are in reasonably 
good agreement. The experiments in 2D are done on 
triangular lattices, and we also study numerically and 
analytically pulse propagation on such lattices with and 
without friction. The experiments for steel spheres and 
the predictions for frictional lattices are in good agree- 
ment (even though there are some subtle points regarding 
the scaling with pressure, see below). 

We also compare our numerical results for the disor- 
dered system with predictions following from numerical 
estimates of the effective long wavelength elastic con- 
stants of our packings. Remarkably, even though elas- 
tic constants predict the time-of-flight velocity reason- 
ably well, there is a 40 % discrepancy between predicted 
and observed velocity for our systems. A possible rea- 
son for this is that our pulses may probe the system on 
short scales which are not governed by a long wavelength 
expansion; indeed a (preliminary) analysis of the spa- 
tial structure of the modes that contribute significantly 
to pulse propagation indicates that most modes appear 
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rather different from simple plane waves. It is likely, but 
very hard to check numerically (at least with the meth- 
ods used in this paper), that for propagation over larger 
distances (such as those probed in the engineering litera- 
ture) the elastic approximation becomes better, and the 
dominant modes would become simple plane waves. The 
crucial open question becomes thus what sets the length 
scale at which such description becomes applicable. Re- 
cently this issue has also emerged within the context of 
the proposal that the static behavior of granular packings 
of hard particles is governed by a critical point ( "point 
J") in a jamming phase diagram Q. We will come back 
to the relation to the jamming phase diagram elsewhere 

We also found that the amplitude of the coherent part 
and its width scales with a power of the distance as the 
signal propagates. For the initial condition where grains 
touching one wall have nonzero velocity, the amplitude 
exponent is 7 « 1.5 for disordered 2D systems, while 
it is 2/3 (exact result) for a ID chain of identical balls. 
The exponent of the signal width is q « 1 for the disor- 
dered 2D system, and 1/3 for the ID chain. The shape 
of the signal can be computed as well, and it is given by 
Airy functions for the ID chain. A triangular lattice of 
identical balls can be mapped (except for the strength of 
the wall springs) to the ID chain, predicting the same 
exponents and signal shape. 

A final issue that we studied in detail is the varia- 
tion of the sound velocity with pressure, since this is an 
important experimental parameter. Our simulations for 
frictional contacts recover the expected p^/^ behavior for 
the time-of-flight velocity and bulk modulus, but not for 
the shear modulus: we found that the transversal wave 
speed scales approximately as p^^"^. These results should 
be compared to results for frictionless sphere packings 
with Hertzian contacts as studied by O'Hcrn et al. 
They found that the bulk modulus B ^ p'^l'^ at low pres- 
sure, while the shear modulus G scaled as p^l^ ^ resulting 
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Some of the experimental data for ctof 
or for resonance frequencies in bead assemblies, and 
some evaluations of elastic moduli in numerical simula- 
tions 0,01, evidence a larger exponent, or at least some 
departure from the p^^^ scaling. The physical origin of 
such observations has been the subject of considerable 
debate In fact, results for the pressure de- 

pendence of sound velocity in disordered glass bead pack- 
ings are somewhat different according to the conditions of 
the experiment, and apparent values of exponents a in a 
c ~ fit vary roughly between 0.16 and 0.25. This calls 
for detailed investigations of the influence of the internal 
state of packings on sound velocities and their pressure 
dependence. While the data published in 0, shown on 
Fig. |S1 indicate a crossover from p^^^ at low p to p^l^ at 
higher pressure, Domenico's results [sof . corresponding 
to much larger confining stress, are fitted by a p^/* law. 



Other data by Jia and Mills (see e.g., ref. 0, figure 
10) agree with c ~ p0-2i ^jjg whole studied pressure 
range, while Sharifipour et al. report in some cases 
exponents a as high as 0.28. Pressure dependences with 
exponents a ~ 0.25 often observed with sands (see e.g., 
[13) are likely to be related to the non-Hertzian behavior 
of contacts between angular particles (or between asper- 
ities of rough particles) as discussed by Goddard j^] , and 
are outside the scope of our simulations. 

Another suggested origin for a different effective scal- 
ing for Hertzian contacts is the increase of coordination 
number with pressure H H H E El, which gradually 
stiffens the packing. In our case, this increase is rather 
small (from z* ~ 3.2 to z* ~ 3.5) and does not entail any 
deviation from the p^^^ scaling for the effective longitu- 
dinal speed Ctof. 

Such an explanation by pressure-induced recruitment 
of additional contacts seems more plausible in regular 
lattices of nominally identical spheres 0, 0, 0, |49j , in 
which a slight polydispersity (or lack of sphericity) causes 
lattice imperfections and strongly reduces the coordina- 
tion number, which only recovers the perfect lattice value 
at high enough confining pressure, P* > 1 (see (|IV Dll ). 
We can expect such a mechanism to explain the experi- 
mental observations by Duffy and Mindlin 49] and Gilles 
and Coste [l3 , as numerical studies of elastic moduli j 
as well as a self-consistent "effective medium" approach 
by Velicky and Caroli 0, both find deviations from a 
p^/^ scaling of long wavelength sound. This effect is of 
course absent in our simulations of perfect lattices. 
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APPENDIX A: CONTACT FORCES IN 
NUMERICAL MODEL 

The normal force between two particles in contact is 
given by the 3D Hertz law [HJl , which is the force between 
two elastic spheres (labeled 1 and 2): 



(Al) 
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where the effective radius R12 = [{Ri) ^ + {R2) ^] ^ 
and effective Young modulus £^1*2 = [(^^i)"^ + [E^y^]-^ 
are half of the harmonic averages of the two grain's pa- 
rameters. Here we introduced the material parameter 
E* = E/{l—v'^) (modified Young modulus, non-standard 
notation), where E is the Young modulus and v is the 
Poisson number. The distance of approach (or "virtual 
normal overlap" ) is given hy n = Ri + R2 — ri2 , where 
ri2 is the distance of the two particle centers. Grain- wall 
interaction can be obtained setting i?wan = 00, and we 
used hard walls (i^waii — 00) ■ 

Implementing the frictional force is less straightfor- 
ward, because frictional contacts can have a memory 
of their history. The standard approach is to consider 
changes in the tangential force with Mindlin's approxi- 
mation I5II: 



(A2) 



'12 

where the elastic constant Gl 
can be calculated from the two grains' material parame- 
ter G* = E/[2{l + v){2- ly)] = E*{1 - v)/{4: - 2v), and 
the virtual tangential displacement t of the particle sur- 
faces is determined from their centers' motion and their 
rotations. This incremental force law is augmented with 
the Coulomb condition: 



|i^t| < ^^Fn^ 



(A3) 



where we take the friction coefficient /i as parameter. It 
is interesting to note that for a given contact the ratio 
of the normal stiffness to the tangential stiffness (despite 
the very different force laws) is constant, we call it 77: 



dFt/dt _ 8G* 
dFn/dn ~ 2E* 



= 1 



(A4) 



assuming the two grains have the same elastic parame- 
ters. For example for v — 0, the value we use in most of 
the simulations, the two effective stiffnesses dF^/dt and 
dFn/dn are equal, so for vibrations neither the radial 
nor the frictional part of the contact will dominate the 
other. (Note that Eq. (|A4|) . a consequence of approxima- 
tion l|A2|l . is really only valid for |-Ft| ^ M-Pn, as pointed 
out in |3(|.) 

The Coulomb condition introduces dissipation, be- 
cause contact surfaces may slip at nonzero force. The 
ensuing dissipation occurs only when the yielding thresh- 
old is exceeded, and not in the infinitesimal amplitude 
oscillations we study here. Nevertheless, in some cases 
we wish to add dissipation, for example when creating 
the packing from the granular gas, or when studying the 
effect of damping on the small amplitude oscillations. 
For this purpose we chose a particular form of damp- 
ing, which is at every instant a given constant fraction 
of the linear critical viscous damping, both for the nor- 
mal and for the tangential force. Through this proce- 
dure, the effective damping force, like the normal and 
tangential Hertz-Mindlin forces, depends nonlinearly on 



the distance of approach n. We impose that the total 
radial force, which now also contains the viscous contri- 
bution, never becomes attractive. 

This choice is appealing theoretically because it con- 
tains only one non-dimensional parameter to control the 
strength of the dissipation. In practice it is not clear 
what the best approximation is of the real (dry or wet 
grains') dissipation. In any case there should be some 
contacts that dominate the dissipation. For those con- 
tacts the viscous force is a certain fraction of the critical 
damping. Our approach is that we impose this ratio on 
all contacts. 



APPENDIX B: SMALL SCALE STATIC 
ELASTICITY 

Let us consider a homogeneous macroscopic sample of 
an elastic material with the same symmetries and bound- 
ary conditions as our numerical systems, and apply a 
body force (per unit volume): 



/(j')(y) = /osin(V^) 



(Bl) 



depending on coordinate y, and directed parallel to axis 
a (1 or 2). With boundary conditions u = on the 
top (y = L2) and bottom (y = 0) walls, and lateral 
periodicity {Li is the system width) the corresponding 
displacement field only has a non-vanishing coordinate 
Ma, given by: 



,(") 



/ X /oi2 ■ /"Try 
" 2 2ry ^^"(^~)' 



with Cq = C22 for a — 2 and Cq 
the total elastic energy 



C33 for a — 1. Hence 



(B2) 



To mimic the force field of Fan. IBll in our discrete sam- 
ples, each non-rattler bead i is submitted to a force : 



AttRI 
3$? 



/o sin( 



m:yi 



while the bottom and top walls are fixed and an apparent 
elastic modulus is obtained from the total elastic energy 
using Fqn. IIj2l The resulting "local" constants, denoted 
as C22{n) and Czzin) are compared in table llVl for n = 1 
and n = 2, to the corresponding "global" values (given 
in table The longitudinal constant C22, measured in 
this way, is only slightly lower, but roughly agrees with 
the previous result. Results for transverse constants are 
similar, except that sample to sample fluctuations are 
somewhat larger. In view of the small sample size, and 
the importance of boundary effects, it can be concluded, 
especially for the higher pressures, that the static elastic 
response to force fields is in reasonable agreement with 
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p 


C22(l)/C22 C22(2)/C22 C33(l)/C33 C33(2)/C33 


10-' 


0.85 ±0.05 0.74 ±0.09 0.81 ±0.14 1.01 ±0.24 


10'^ 


0.88 ± 0.04 0.80 ± 0.06 0.90 ± 0.20 0.96 ±0.16 


10"^ 


0.91 ± 0.05 0.86 ± 0.06 0.85 ±0.10 0.89 ±0.11 


10-* 


0.93 ±0.04 0.88 ±0.04 0.85 ±0.10 0.90 ± 0.07 



TABLE IV: Ratio of apparent elastic moduli deduced from 
1B2I I for n — 1 and n = 2 to the "global" values. 



the equations of elasticity involving the moduli measured 
on globally deforming the sample. This is further con- 
firmed by another set of static response calculations, in 
which displacements, rather than forces are imposed. Let 
us define regularly spaced horizontal lines through the 
sample at y = kL2/(n + 1) for < fc < n ± 1 (so that 
k = Q corresponds to the bottom and k — n + \ corre- 
sponds to the top), with n an odd number. Let us impose 
constant displacements u = in direction a on line k if 
k is even, u = (— l)'uo on line fc if fc = 2Z ± 1 is odd. 
In a homogeneous elastic system, the displacement field 
varies linearly between neighboring lines k and /e + 1, and 
the total elastic energy reads 



(B3) 



Imposing to the center of each particle crossed by the k 
lines the displacements of the same point in a homoge- 
neous continuum (and leaving the rotations free), com- 
puting the total elastic energy and using formula lB3l one 
deduces other values for C22 and C33, denoted as (^^(n), 
given in table for n — \ and n = 3. Once again, 
they are fairly close to the "global" moduli of table ITl|l . 
especially for the higher p values, albeit slightly larger. 
We therefore conclude that the elastic moduli are only 



V 


C'22(l)/C22 C22(3)/C22 C'33(l)/C33 C'33(3)/C33 


10-^ 


1.01 ±0.02 1.09 ± 0.03 1.26 ±0.16 2.49 ± 0.69 


10-*^ 


1.01 ±0.01 1.07 ±0.07 1.28 ±0.20 2.05 ± 0.45 


10-^ 


1.01 ±0.02 1.06 ±0.03 1.15 ±0.11 1.65 ±0.30 


10-* 


1.00 ±0.01 1.03 ±0.02 1.08 ±0.05 1.37 ±0.14 



TABLE V: Ratio of apparent elastic moduli deduced from 
1B3II for n = 1 and n = 3 to the "global" values. 

weakly affected by finite size effects. 



APPENDIX C: ID CHAIN OF IDENTICAL 
BALLS 

We model the ID chain of identical balls with TV identi- 
cal particles of mass to, separated by distance and con- 
nected by linear springs of stiffness K. The first and last 
ball is also connected with identical springs to walls. This 
system models small amplitude oscillations of Hertzian 



balls under finite static pressure. The n-th eigenmode of 
this simple linear system is given by 



u^x^{t) — sin(fc„a;) sin(a;„t + 0„) . 



(CI) 



where u{t) is the displacement of a ball, and we label the 
balls by their position x = i£, where i ~ 1,2, . . . , N. The 
wavenumbers and eigenfrequencies are 

1 '^'^ ^ [K . kni 

The above dispersion relation determines the longitudi- 
nal sound speed 



duj 
dk 



k=0 



(C3) 



by which the long wavelength waves propagate. The full 
solution is given by 



N 



"^xit) = ^ a„ sin(fc„x) sin(w„t ± 0„) 



(C4) 



where the amplitudes a„ are computed by projecting the 
initial condition onto the modes. The two cases consid- 
ered here are (A) ui{t = 0) = c, and (B) ui{t = Q) = t, 
all other displacements and velocities at i = are zero. 
This gives 



N 



2^ kfi^ 

cos 

1 2 



2£ 



iV± 1 



sinfc„£. (C5) 



We are interested in the time dependence of the force the 
A^-th ball presses the wall: 



F{t) = KuN{t) 



(C6) 



We cannot calculate this in closed form, but can for ex- 
ample look at its value at the time "a wave would arrive" : 
to = Nljc. Substituting the first initial condition's 
into Eq. HC4() and rewriting the highly oscillating terms 
we get 



N 



2Ki kni . , 

— COS — — sm kn 

N + 1 2 

X sin knl + 

V ^0 



(C7) 



o{km 



1,3 



where fcg = (24/iV)^/'^/£. The terms in the sum become 
highly oscillating for fc„ > k^, effectively canceling each 
other. The dominant contribution therefore comes from 
terms with fc„ < fcg (or equivalently kn < constant x feg, 
as ultimately we will only compute scaling exponents). 
The sum of the slowly varying terms is approximated by 
integral: 



ko 



2Kedk ke 
cos — sm k£ 

TT 2 



(C8) 
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We look at the asymptotics N ^ oo (implying fco — > 0), The above two relations immediately lead to Eqs. (j ^lO|l . 
for which we have to find the terms lowest order in fep- 
This yields 

F^ito) ~ kl (C9) 

Interestingly the second initial condition gives a different 
answer: 

i^^(io) ~ kl (CIO) 
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